#include "dEdxGausN.h"

void fith()
{
  TList *ll = new TList;
  Int_t col[]={kBlue, kRed, kGreen+3, kMagenta, kYellow};
    
  //________________________________________________________________
  //                   initialize input
  //________________________________________________________________
  
  //original without any treatment
  TH2D *hkk = (TH2D*)gDirectory->Get("kk");// pt:tpcip
  if(!hkk){
    printf("fith error no kk!\n"); gDirectory->ls(); exit(1);
  }
  ll->Add(hkk);

  TH2D *htoflogp = (TH2D*)gDirectory->Get(Form("h%slogp", dEdxGausN::TOFTag()));
  if(!htoflogp){
    printf("fith error no toflogp. in noTOF mode!\n"); gDirectory->ls();
  }
  else{
    ll->Add(htoflogp);
  }

  const TString horiName(Form("h%soriginal", dEdxGausN::TOFTag()));
  TH2D *horiginal =  (TH2D*)gDirectory->Get(horiName);
  printf("fith horiginal name: %s\n", horiName.Data());
  if(!horiginal){
    printf ("no horiginal!\n"); gDirectory->ls(); exit(1);
  }
  ll->Add(horiginal);
  
  //original hyield without any treatment
  TH1D *horiginalyield =  (TH1D*)gDirectory->Get(horiName+"yield");
  if(!horiginalyield){
    printf ("no horiginalyield!\n"); gDirectory->ls(); exit(1);
  }
  ll->Add(horiginalyield);

  TGraphErrors * greftot=(TGraphErrors*)gDirectory->Get("greftot");
  if(dEdxGausN::IsWithRef() && !greftot){
    printf("fith dEdxGausN::IsWithRef() && !greftot %d\n", dEdxGausN::IsWithRef()); gDirectory->ls(); exit(1);
  }
  if(greftot){
    ll->Add(greftot);
  }

  TGraphErrors * gref[dEdxGausN::NparType()];
  for(Int_t itype=0; itype<dEdxGausN::NparType(); itype++){
    gref[itype]=(TGraphErrors*)gDirectory->Get(Form("gref%d", itype));
    if(gref[itype]){
      ll->Add(gref[itype]);
    }
  }

  TH2D *hfitdata = 0x0;
  TH1D *hdummyyield = 0x0;
  dEdxGausN::GetHFit(horiginal, horiginalyield, hfitdata, hdummyyield);
  ll->Add(hfitdata);
  ll->Add(hdummyyield);

  //_________________________________________________________________
  //                      initialize numbers and controls
  //_________________________________________________________________

  Double_t bbpar[100];
  memcpy(bbpar, dEdxGausN::IniBBpar(), sizeof(Double_t)* dEdxGausN::NparBB());
  for(Int_t ii=0; ii<dEdxGausN::NparBB(); ii++){
    printf("IniBBpar %d: %f\n", ii, bbpar[ii]);
  }

  const Int_t npar = dEdxGausN::Npar();
  Double_t pars[npar], errs[npar], chi[100], lows[npar], highs[npar];
  chi[0]=chi[1]=-999;
  Bool_t pfix[npar];
  Int_t order[dEdxGausN::Ntype()];
  
  for(Int_t ii=0; ii<npar; ii++){
    errs[ii] = -999;
    lows[ii]= 0;
    highs[ii]= 0;
  }
  
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    //ini value
    pars[itype*dEdxGausN::NparType()+0]= 0.5;
    pars[itype*dEdxGausN::NparType()+1]= -999;
    pars[itype*dEdxGausN::NparType()+2]= dEdxGausN::IniResolution();

    //limits
    lows[ itype*dEdxGausN::NparType()+0]  = 0;
    highs[itype*dEdxGausN::NparType()+0]  = 1;

    lows[ itype*dEdxGausN::NparType()+2]  = dEdxGausN::MinResolution();
    highs[itype*dEdxGausN::NparType()+2]  = dEdxGausN::MaxResolution();

    //order
    order[itype] = itype;
  }

  //is fixed
  dEdxGausN::ResetIfixed(pfix);

  //________________________________________________________________
  //              check old aleph fit from directory
  //________________________________________________________________
  TList *lgr = new TList;

  Int_t irev = 0;
  TH1D * oldmean = 0x0;

  TList *tmplg = (TList*) gDirectory->Get("lgr");

  while(tmplg){
    TH1D * tmpom = (TH1D*) tmplg->FindObject(Form("alephFITRev%d", irev));

    if(tmpom){
      lgr->Add(tmpom);
      oldmean = tmpom;

      TString grnames[100]={"allp1VSbg"};
      Int_t ngr=1;
      for(Int_t tmpit=0; tmpit<dEdxGausN::Ntype(); tmpit++){
        grnames[ngr++]=Form("allp1VSmomType%d", tmpit);
        grnames[ngr++]=Form("allp2VSp1Type%d", tmpit);
        grnames[ngr++]=Form("allp2VSmomType%d", tmpit);
      }

      for(Int_t igr=0; igr<ngr; igr++){
        const TString grn(Form("%sRev%d",grnames[igr].Data(), irev));
        TGraphErrors * tmpgr = (TGraphErrors*) tmplg->FindObject(grn);
        if(!tmpgr){
          printf("fith no %s %d while tmpom exist!\n", grnames[igr].Data(), irev); tmplg->ls(); exit(1);
        }
        lgr->Add(tmpgr);
      }
      
      irev++;
    }
    else{
      break;
    }
  }

  printf("fith rev %d not found, fit rev %d now!\n", irev, irev);

  if(oldmean){
    TString omtit(oldmean->GetTitle());
    for(Int_t iom=0; iom<dEdxGausN::NparBB(); iom++){
      const TString tmpf = omtit(0,omtit.First(',')-1);
      bbpar[iom]=tmpf.Atof();
      printf("fith old version mean used %s: %d %e\n", oldmean->GetTitle(), iom, bbpar[iom]);
      omtit=omtit(omtit.First(',')+1,omtit.Length());
    }
  }

  //________________________________________________________________
  //           initialize output, not all marked with irev
  //________________________________________________________________
  
  TGraphErrors * allp1VSbg = new TGraphErrors;
  allp1VSbg->SetName(Form("allp1VSbgRev%d",irev));
  lgr->Add(allp1VSbg);
  
  TGraphErrors * allp1VSmomTypes[dEdxGausN::Ntype()];
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    allp1VSmomTypes[itype] = new TGraphErrors;
    allp1VSmomTypes[itype]->SetMarkerColor(col[itype]);
    allp1VSmomTypes[itype]->SetMarkerSize(1);
    allp1VSmomTypes[itype]->SetMarkerStyle(20+itype);
    allp1VSmomTypes[itype]->SetName(Form("allp1VSmomType%dRev%d",itype, irev));
    lgr->Add(allp1VSmomTypes[itype]);
  }

  TGraphErrors * allp2VSp1Types[dEdxGausN::Ntype()];
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    allp2VSp1Types[itype] = new TGraphErrors;
    allp2VSp1Types[itype]->SetMarkerColor(col[itype]);
    allp2VSp1Types[itype]->SetMarkerSize(1);
    allp2VSp1Types[itype]->SetMarkerStyle(20+itype);
    allp2VSp1Types[itype]->SetName(Form("allp2VSp1Type%dRev%d",itype, irev));
    lgr->Add(allp2VSp1Types[itype]);
  }

  TGraphErrors * allp2VSmomTypes[dEdxGausN::Ntype()];
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    allp2VSmomTypes[itype] = new TGraphErrors;
    allp2VSmomTypes[itype]->SetMarkerColor(col[itype]);
    allp2VSmomTypes[itype]->SetMarkerSize(1);
    allp2VSmomTypes[itype]->SetMarkerStyle(20+itype);
    allp2VSmomTypes[itype]->SetName(Form("allp2VSmomType%dRev%d",itype, irev));
    lgr->Add(allp2VSmomTypes[itype]);
  }

  TTreeSRedirector *fStream = new TTreeSRedirector(Form("outfith_%sRev%d.root", dEdxGausN::Opt(), irev));
  
  TList *lfit = new TList;

  const TString treename(Form("fith"));
  //________________________________________________________________
  //________________________________________________________________
  
  Int_t ifit = 0;
  Int_t nfail = 0;

  for(Int_t ix = 1; ix<=hfitdata->GetNbinsX(); ix++){
    Double_t xx = -999;
    xx = hfitdata->GetXaxis()->GetBinCenter(ix);
    if(xx<dEdxGausN::Xmin())
      continue;

    const Double_t itot = horiginalyield->GetBinContent(ix);
    printf("ix %d xx %e ifit %d tot %e\n", ix, xx, ifit, itot);

    if(ix<dEdxGausN::StartBin() || itot<=dEdxGausN::PbinThres()){
      continue;
    }

    if(dEdxGausN::FilterIX(ix))
      continue;

    TH1D *tmph = hfitdata->ProjectionY(Form("%sx%d", hfitdata->GetName(), ix),ix,ix,"oe");
    tmph->SetTitle(Form("x %e",xx));
    lfit->Add(tmph);

    //==============================================================================
    //  fit starts ...
    //==============================================================================
    //oder[ii] is the species in order ii
    Int_t oldid[dEdxGausN::Ntype()];
    memcpy(oldid, order, sizeof(oldid));
    dEdxGausN::SetMean(xx, bbpar, pars, order);
    for(Int_t ipartype=0; ipartype<dEdxGausN::NparType(); ipartype++){
      if(ipartype==1){
        continue;
      }
      dEdxGausN::Swap(xx, ipartype, pars, oldid, order);
    }

    Int_t kfail = BaseUtils::BinnedLikelihoodFit(tmph, dEdxGausN::Func, npar, pars, errs, pfix, 0x0, lows, highs);

    //the mean of the first peak should not collaspe
    if(pars[1]<tmph->GetBinLowEdge(2))
      kfail = 2;

    if(kfail)
      nfail++;

    printf("fith fitresult irev %d ix %d kfail %d chi2ndof %.1f %.0f %.2f\n", irev, ix, kfail, chi[0], chi[1], chi[1]>EPSILON?chi[0]/chi[1]:-999);
    //==============================================================================
    //  fit finished.
    //==============================================================================

    Double_t totalyield = horiginalyield->GetBinContent(ix);
    //tmph->Scale(totalyield);
    tmph->SetLineColor(kBlack);
    //printf("longtestpars_0 %d %e %e %e %e\n",ix, pars[0], pars[3], errs[0], errs[3]);

    Double_t yields[dEdxGausN::Ntype()];
    Double_t outpar[npar];    
    Double_t outerr[npar];
    Int_t ordinv[dEdxGausN::Ntype()];
    dEdxGausN::SavePar(fStream, treename, &totalyield, yields, pars, errs, errs, order, outpar, outerr, outerr, ordinv);
    //printf("longtestpars_1 %d %e %e %e %e\n",ix, pars[0], pars[3], errs[0], errs[3]);

    //fill allp1VSbg, allp1VSmomTypes[], allp2VSp1Types[], allp2VSmomTypes[]
    Double_t bg[dEdxGausN::Ntype()];
    //bg$(itype) should appear in the same order for all ix, so only itype as iter possible
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      bg[itype] = TMath::Power(10, xx -TMath::Log10(dEdxGausN::Mass(itype)));
      (*fStream)<<treename<<
        Form("bg%d=", itype)<<bg[itype];

      //electron dEdx
      //10h 75~90
      //h0e 75~80
      //bg<5e3 to avoid the lowering due to dedx crossing
      const Int_t ith = ordinv[itype];
      const Double_t tmpp1 = outpar[ith*dEdxGausN::NparType()+1];
      const Double_t tmpp2 = outpar[ith*dEdxGausN::NparType()+2];
      //have to be suitablly large to ensure successful of dedxfit, so that the initial bbpar will be updated physically!!
      //in later iteration, realistic error can be used
      //const Double_t tmpe1 = (irev>=dEdxGausN::LastfithIter()-1 && (bg[itype]<1.5 || bg[itype]>200) )? outerr[ith*dEdxGausN::NparType()+1] : 500/pow(10,irev);
      const Double_t tmpe1 = 500/pow(10,irev);
      const Double_t tmpe2 = 1e-4;

      const Int_t ip1VSmom = allp1VSmomTypes[itype]->GetN();
      allp1VSmomTypes[itype]->SetPoint(ip1VSmom, xx, tmpp1);
      allp1VSmomTypes[itype]->SetPointError(ip1VSmom, 0, tmpe1);
      
      //const Double_t tmpmom = TMath::Power(10, xx);
      if( kfail==0 && ( (itype!=dEdxGausN::kELECTRON && tmpp1 < 140) || (bg[itype]<5e3 && tmpp1 > 75 && tmpp1 < 90) ) ){
        const Int_t ip1VSbg = allp1VSbg->GetN();
        allp1VSbg->SetPoint(ip1VSbg, bg[itype], tmpp1);
        allp1VSbg->SetPointError(ip1VSbg, 0, tmpe1);
        
        const Int_t ip2VSp1 = allp2VSp1Types[itype]->GetN();
        allp2VSp1Types[itype]->SetPoint(ip2VSp1, tmpp1, tmpp2);
        allp2VSp1Types[itype]->SetPointError(ip2VSp1, tmpe1, tmpe2);
        
        //see curvature dependence of resolution at small range of mean
        if(tmpp1 > 40 && tmpp1 < 90){//with 40-70 basically only proton, kaon and pion are seen
          const Int_t ip2VSmom = allp2VSmomTypes[itype]->GetN();
          allp2VSmomTypes[itype]->SetPoint(ip2VSmom, xx, tmpp2);
          allp2VSmomTypes[itype]->SetPointError(ip2VSmom, 0, tmpe2);       
        }
      }
    }
    
    Double_t chidof = chi[0]/chi[1];
    (*fStream)<<treename<<
      "irev="<<irev<<
      "ix="<<ix<<
      "xx="<<xx<<
      "chi2="<<chi[0]<<
      "ndof="<<chi[1]<<
      "chidof="<<chidof<<
      "kfail="<<kfail<<
      "ifit="<<ifit<<
      "\n";

    TH1D *hfit = BaseUtils::GetHfit(Form("%sfit", tmph->GetName()), dEdxGausN::Func, pars, tmph->GetXaxis()->GetXmin(), tmph->GetXaxis()->GetXmax());
    hfit->SetLineColor(kBlack);
    //hfit->Scale(totalyield);
    //hfit->SetDrawOption("same hist C");
    lfit->Add(hfit);
    
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      dEdxGausN::DrawTF1(itype, tmph->GetName(), &(outpar[ordinv[itype]*dEdxGausN::NparType()]), col[itype], lfit); 
    }

    //printf("longtestpars_2 %d %e %e %e %e\n",ix, pars[0], pars[3], errs[0], errs[3]);
    //printf("longtestpars_3\n");

    ifit++;
  }
  
  if(nfail>ifit/2){
    printf("fith nfail>ifit/2!! indicating very bad ini bbpar!! %d %d \n", nfail, ifit);exit(1);
  }

  //fit gallmean
  BaseUtils::SetFitPrintLevel(0);
  //if using ref (for irev>=1), then all versions should be identical
  const Bool_t kfitref = dEdxGausN::IsWithRef() && irev;
  //irev0  should be  successful fit in order to update ini bbpar!! -- not really
  Double_t fithlows[10]={0,0,0,0,0,0,0,0,0,0};
  Double_t fithhighs[10]={0,0,0,0,0,0,0,0,0,0};
  if(dEdxGausN::BBOpt()==dEdxGausN::kRawALEPH){
    //fithlows[0]  = 78;
    //fithhighs[0] = 82;
    //fithlows[3]  = 2.3;
    //fithhighs[3] = 2.5;
  }
  const Double_t chi2max = 100;
  TH1D *fitallmean = dEdxGausN::dEdxFit( kfitref ? greftot : allp1VSbg, Form("alephFITRev%d", irev), kFALSE, bbpar, fithlows, fithhighs, chi2max);
  lgr->Add(fitallmean);

  fStream->GetFile()->cd();
  ll->Write();
  lfit->Write("lfit",TList::kSingleKey);
  lgr->Write("lgr",TList::kSingleKey);

  delete fStream;
}

int main(int argc, char* argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2\n"); exit(1);
  }

  const TString fin(argv[1]);
  TFile::Open(fin);

  dEdxGausN::Ini();
  
  fith();

  return 0;
}
